{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Manybody Potential Expressions\n",
    "\n",
    "This notebook uses symbolic differentiation to derive expressions for common manybody potentials defined in terms of the following functions:\n",
    "\n",
    "$$ \\phi(r_{ij}^2, \\xi_{ij}), \\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2) $$\n",
    "\n",
    "To avoid verbosity, expressions in this notebook will use the following symbol substitution:\n",
    "\n",
    "$$ R = R_1 = r_{ij}^2\\\\\n",
    "R_2 = r_{ik}^2\\\\\n",
    "R_3 = r_{jk}^2\\\\\n",
    "\\xi = \\xi_{ij} $$\n",
    "\n",
    "So that the above functions become:\n",
    "\n",
    "$$ \\phi(R, \\xi), \\Theta(R_1, R_2, R_3) $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from sympy.abc import R, xi, x\n",
    "from IPython.display import display, Latex\n",
    "\n",
    "init_printing()\n",
    "\n",
    "R1, R2, R3 = symbols(\"R_{1:4}\")\n",
    "a, x = symbols(\"a, x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are two classes that represent the above functions $\\phi$ and $\\Theta$. Their member functions `gradient()` and `hessian()` return the correct derivatives in the order expected by the `Manybody` calculator class in `matscipy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Potential:\n",
    "    def __init__(self, expression, latex):\n",
    "        self.e = expression\n",
    "        self.latex = latex\n",
    "    def show(self):\n",
    "        display(Latex(f\"$${self.latex} = {latex(self.e)}$$\"))\n",
    "        display(Latex(f\"$$\\\\nabla{self.latex} = {latex(self.gradient)}$$\"))\n",
    "        display(Latex(f\"$$\\\\nabla^2{self.latex} = {latex(self.hessian)}$$\"))\n",
    "\n",
    "class Phi(Potential):\n",
    "    def __init__(self, expression):\n",
    "        super().__init__(expression, \"\\\\phi(R, \\\\xi)\")\n",
    "        self.gradient = [(diff(self.e, v)) for v in (R, xi)]\n",
    "        self.hessian = [(diff(self.e, *v)) for v in [(R, R), (xi, xi), (R, xi)]]\n",
    "\n",
    "class Theta(Potential):\n",
    "    def __init__(self, expression):\n",
    "        super().__init__(expression, \"\\\\Theta(R_1, R_2, R_3)\")\n",
    "        self.gradient = [(diff(self.e, v)) for v in (R1, R2, R3)]\n",
    "        self.hessian = [(diff(self.e, *v)) for v in [(R1, R1), (R2, R2), (R3, R3), (R2, R3), (R1, R3), (R1, R2)]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Harmonic Pair Potential\n",
    "\n",
    "Below is the definition of a harmonic pair potential, with the added contribution from the three body interaction $\\xi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\phi(R, \\xi) = \\frac{k \\left(\\sqrt{R} - r_{0}\\right)^{2}}{2} + \\xi$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla\\phi(R, \\xi) = \\left[ \\frac{k \\left(\\sqrt{R} - r_{0}\\right)}{2 \\sqrt{R}}, \\  1\\right]$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla^2\\phi(R, \\xi) = \\left[ \\frac{k \\left(\\frac{1}{R} - \\frac{\\sqrt{R} - r_{0}}{R^{\\frac{3}{2}}}\\right)}{4}, \\  0, \\  0\\right]$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define equilibrium distance and stiffness\n",
    "k, r0 = symbols(\"k, r_0\")\n",
    "\n",
    "energy = k / 2 * (sqrt(R) - r0)**2 + xi\n",
    "harmonic = Phi(energy)\n",
    "harmonic.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lennard-Jones Pair Potential\n",
    "\n",
    "LJ works nicely because of even powers that play nice with the definition in terms of the *squared distance*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\phi(R, \\xi) = 4 \\epsilon \\left(- \\frac{\\sigma^{6}}{R^{3}} + \\frac{\\sigma^{12}}{R^{6}}\\right) + \\xi$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla\\phi(R, \\xi) = \\left[ 4 \\epsilon \\left(\\frac{3 \\sigma^{6}}{R^{4}} - \\frac{6 \\sigma^{12}}{R^{7}}\\right), \\  1\\right]$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla^2\\phi(R, \\xi) = \\left[ - \\frac{24 \\epsilon \\sigma^{6} \\cdot \\left(2 - \\frac{7 \\sigma^{6}}{R^{3}}\\right)}{R^{5}}, \\  0, \\  0\\right]$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define characteristic distance and energy\n",
    "eps, sigma = symbols(\"\\\\epsilon, \\\\sigma\")\n",
    "\n",
    "energy = 4 * eps * (sigma**12 / R**6 - sigma**6 / R**3) + xi\n",
    "lj = Phi(energy)\n",
    "lj.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Harmonic Angle Potential\n",
    "\n",
    "Three-body potential that enforces an angle stiffness between two bonds.\n",
    "\n",
    "$$ \\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2) = \\frac{1}{2} k_\\Theta (\\Theta_{ijk} - \\Theta_0)^2 = E(a)$$\n",
    "\n",
    "with \n",
    "\n",
    "$$ E =  \\frac{1}{2} k_\\Theta (a - \\Theta_0)^2$$\n",
    "$$ a = \\arccos(x)$$\n",
    "\n",
    "The gradient with respect to distance/ squared distance $r_x$/$R_x=r_x^2$ with $x=\\{ij, ik, jk\\}$ reads\n",
    "\n",
    "$$\\frac{\\partial\\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2)}{\\partial r_x} = \\frac{\\partial \\Theta}{\\partial x} \\frac{\\partial x}{\\partial r_x}$$\n",
    "\n",
    "$$\\frac{\\partial\\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2)}{\\partial R_x} = \\frac{\\partial \\Theta}{\\partial x} \\frac{\\partial x}{\\partial r_x} \\frac{\\partial r_x}{\\partial R_x}$$\n",
    "\n",
    "The Hessian with respect to distance/ squared distance reads\n",
    "\n",
    "$$\n",
    "\\frac{\\partial^2 \\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2)}{\\partial r_y \\partial r_x} \n",
    "= \n",
    "\\frac{\\partial^2 \\Theta}{\\partial y \\partial x}\\frac{\\partial y}{\\partial r_y} \\frac{\\partial x}{\\partial r_x}\n",
    "+\n",
    "\\frac{\\partial \\Theta}{\\partial r_x} \\frac{\\partial^2 x}{\\partial r_y \\partial r_x}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial^2 \\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2)}{\\partial R_y \\partial R_x} \n",
    "= \n",
    "\\left(\n",
    "\\frac{\\partial^2 \\Theta}{\\partial y \\partial x}\\frac{\\partial y}{\\partial r_y} \\frac{\\partial x}{\\partial r_x}\n",
    "+\n",
    "\\frac{\\partial \\Theta}{\\partial r_x} \\frac{\\partial^2 x}{\\partial r_y \\partial r_x}\n",
    "\\right)\n",
    "\\frac{\\partial r_x}{\\partial R_x}\\frac{\\partial r_y}{\\partial R_y}\n",
    "+\n",
    "\\frac{\\partial \\Theta}{\\partial x}\\frac{\\partial x}{\\partial r_x} \\frac{\\partial^2 r_x}{\\partial R_y \\partial R_x}\n",
    "$$ \n",
    "\n",
    "Additional we have the derivatives\n",
    "$$\n",
    "\\frac{\\partial \\Theta}{\\partial x} = \\frac{\\partial \\Theta}{\\partial a} \\frac{\\partial a}{\\partial x}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial^2 \\Theta}{\\partial y \\partial x} \n",
    "=\n",
    "\\frac{\\partial^2 \\Theta}{\\partial a^2} \\frac{\\partial a}{\\partial y} \\frac{\\partial a}{\\partial x}\n",
    "+\n",
    "\\frac{\\partial \\Theta}{\\partial a} \\frac{\\partial^2 a}{\\partial y \\partial x}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\Theta(R_1, R_2, R_3) = \\frac{R_{1} - R_{2} + R_{3}}{2 \\sqrt{R_{1} R_{3}}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla\\Theta(R_1, R_2, R_3) = \\left[ \\frac{1}{2 \\sqrt{R_{1} R_{3}}} - \\frac{R_{1} - R_{2} + R_{3}}{4 R_{1} \\sqrt{R_{1} R_{3}}}, \\  - \\frac{1}{2 \\sqrt{R_{1} R_{3}}}, \\  \\frac{1}{2 \\sqrt{R_{1} R_{3}}} - \\frac{R_{1} - R_{2} + R_{3}}{4 R_{3} \\sqrt{R_{1} R_{3}}}\\right]$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla^2\\Theta(R_1, R_2, R_3) = \\left[ \\frac{-4 + \\frac{3 \\left(R_{1} - R_{2} + R_{3}\\right)}{R_{1}}}{8 R_{1} \\sqrt{R_{1} R_{3}}}, \\  0, \\  \\frac{-4 + \\frac{3 \\left(R_{1} - R_{2} + R_{3}\\right)}{R_{3}}}{8 R_{3} \\sqrt{R_{1} R_{3}}}, \\  \\frac{1}{4 R_{3} \\sqrt{R_{1} R_{3}}}, \\  \\frac{- \\frac{2}{R_{3}} - \\frac{2}{R_{1}} + \\frac{R_{1} - R_{2} + R_{3}}{R_{1} R_{3}}}{8 \\sqrt{R_{1} R_{3}}}, \\  \\frac{1}{4 R_{1} \\sqrt{R_{1} R_{3}}}\\right]$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\Theta(R_1, R_2, R_3) = 0.5 k_{\\theta} \\left(- \\theta_{0} + a\\right)^{2}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla\\Theta(R_1, R_2, R_3) = 0.5 k_{\\theta} \\left(- 2 \\theta_{0} + 2 a\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla^2\\Theta(R_1, R_2, R_3) = 1.0 k_{\\theta}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\Theta(R_1, R_2, R_3) = \\operatorname{acos}{\\left(x \\right)}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla\\Theta(R_1, R_2, R_3) = - \\frac{1}{\\sqrt{1 - x^{2}}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\nabla^2\\Theta(R_1, R_2, R_3) = - \\frac{x}{\\left(1 - x^{2}\\right)^{\\frac{3}{2}}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "k_theta, theta0 = symbols(\"k_\\\\theta, \\\\theta_0\")\n",
    "\n",
    "cos_angle = (R1 + R3 - R2)/(2 * sqrt(R1 * R3))\n",
    "\n",
    "# Print derivatives of angle cosine\n",
    "Theta(cos_angle).show()\n",
    "\n",
    "# Derivatives with respect to x \n",
    "E = 0.5 * k_theta * (a - theta0)**2\n",
    "dTheta_da = Theta(E)\n",
    "dTheta_da.gradient = diff(E, a)\n",
    "dTheta_da.hessian = diff(E, *(a, a))\n",
    "dTheta_da.show()\n",
    "\n",
    "a = acos(x)\n",
    "da_dx = Theta(a)\n",
    "da_dx.gradient = diff(a, x)\n",
    "da_dx.hessian = diff(a, *(x, x))\n",
    "da_dx.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stillinger-Weber Potential\n",
    "\n",
    "The functional form of the Stillinger-Weber potential reads\n",
    "\n",
    "$$ \\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2) = \\frac{1}{2} k_\\Theta (\\Theta_{ijk} - \\Theta_0)^2 = E(a)$$\n",
    "\n",
    "with \n",
    "\n",
    "$$ E =  \\frac{1}{2} k_\\Theta (a - \\Theta_0)^2$$\n",
    "$$ a = \\arccos(x)$$\n",
    "\n",
    "The gradient with respect to distance/ squared distance $r_x$/$R_x=r_x^2$ with $x=\\{ij, ik, jk\\}$ reads\n",
    "\n",
    "$$\\frac{\\partial\\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2)}{\\partial r_x} = \\frac{\\partial \\Theta}{\\partial x} \\frac{\\partial x}{\\partial r_x}$$\n",
    "\n",
    "$$\\frac{\\partial\\Theta(r_{ij}^2, r_{ik}^2, r_{jk}^2)}{\\partial R_x} = \\frac{\\partial \\Theta}{\\partial x} \\frac{\\partial x}{\\partial r_x} \\frac{\\partial r_x}{\\partial R_x}$$\n",
    "\n",
    "The Hessian with respect to distance/ squared distance reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
